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Abstract 

The human intestinal microbiota performs many essential functions for the host. Antimicrobial agents, such as antibiotics 
(AB), are also known to disturb microbial community equilibrium, thereby having an impact on human physiology. While an 
increasing number of studies investigate the effects of AB usage on changes in human gut microbiota biodiversity, its 
functional effects are still poorly understood. We performed a follow-up study to explore the effect of ABs with different 
modes of action on human gut microbiota composition and function. Four individuals were treated with different 
antibiotics and samples were taken before, during and after the AB course for all of them. Changes in the total and in the 
active (growing) microbiota as well as the functional changes were addressed by 16S rRNA gene and metagenomic 454- 
based pyrosequencing approaches. We have found that the class of antibiotic, particularly its antimicrobial effect and 
mode of action, played an important role in modulating the gut microbiota composition and function. Furthermore, analysis 
of the resistome suggested that oscillatory dynamics are not only due to antibiotic-target resistance, but also to fluctuations 
in the surviving bacterial community. Our results indicated that the effect of AB on the human gut microbiota relates to the 
interaction of several factors, principally the properties of the antimicrobial agent, and the structure, functions and 
resistance genes of the microbial community. 
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Introduction 

Throughout evolution mammals have established symbioses 
with microbial communities, which are located in different organs 
and tissues of the body such as skin, mucosa, or the gastrointestinal 
tract. The gut microbiota in humans is a particularly complex 
ecosystem with few dominant phyla (Firmicutes, Bacteroidetes, 
Proteobacteria and Actinobacteria) but show greater microbial 
diversity at lower taxonomic levels and a high functional 
redundancy [1,2]. The gut microbiota seems to be host-specific 
and rather stable under non- or small perturbations [3] and is 
involved in a large number of host beneficial functions such as 
food processing, growth regulation of the intestinal epithelium, 
development of the immune system, or protection against 
pathogens [2,4,5]. Because of the essential role of the microbiota 
in host life, imbalances in the gut microbial community may have 
an important impact on human health. This is apparent in some 
intestinal pathologies such as inflammatory bowel diseases or 
antibiotic-associated diarrhea [6]. 

Systematic antibiotic (AB) therapy represents a major public 
health problem because gut microbiota may be transformed into a 
reservoir of antibiotic resistance genes, promoting the appearance 



of harmful resistant strains [7,8,9,10]. It also suppresses some 
protective members of the resident microbiota promoting over- 
growth of opportunistic pathogens such as Cbstiidium difficile [11]. 
Moreover, AB therapy disturbs the gut microbiota and, concom- 
itantly, affects human physiology, for instance carbohydrate 
metabolism or immunity [7,12]. 

Antibiotic features such as class, spectrum or pharmacological 
properties affect the gut microbiota in different ways [13]. In 
addition, host-associated factors such as diet, life history, genetic or 
health status, properties of the gut microbial ecosystem itself like 
resistance and resilience, or even the interplay between the 
microbiota and its host also have an effect on microbiota 
composition and function. All these factors can mask changes 
caused exclusively by antibiotics, representing a real challenge 
when it comes to understand microbiota responses. Most recent 
studies into the impact of antibiotics on the microbiota have 
focused on the emergence of resistant strains, but few have 
described their influence on the microbial community itself 
[14,15,16,17,18]. These latter surveys, using 16S rRNA gene 
sequencing, have shown that short and long-term AB courses 
affect diversity and biomass of the intestinal microbiota, with 
microbial composition resilience remaining deficient for long time 
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after AB-treatment [15,16,17,18]. By contrast, the functional 
impact of AB on the microbial ecosystem has been addressed less 
frequently [19]. 

The use of the meta-"omics" approaches (metagenomics, 
metatranscriptomics, metaproteomics) has provided deeper in- 
sights into microbial communities in different ecosystems 
[2,20,21,22,23]. A recent integrated analysis has provided a better 
understanding of the nature of the complex processes underlying 
the whole human gut microbiota and its responses during beta- 
lactamic-therapy [12]. 

In the present work we studied the effect of different antibiotics 
on the human gut microbiota by a follow-up study comparing 
microbial communities before, during and after AB therapy in four 
individuals. We analyzed the changes in composition of the total 
(16S rRNA gene) and active (16S rRNA transcripts) microbiota 
throughout treatment. Furthermore, the functional analysis of the 
total gene content of the community showed, for the first time, 
how the mode of action and the antimicrobial effect of AB affected 
the functional potential of the community. Finally, we described 
the dynamics of resistance genes (i.e. the resistome) throughout the 
study, paying particular attention to those that become resident 
after AB-therapy. 

Materials and Methods 

Ethics statement 

The study was approved by the Ethics Committee of the 
Medical Faculty of the Christian-Albrechts University Kiel, 
Germany. Informed written consent was obtained from all patients 
involved in the study. 

Sample collection and AB treatment regimen 

Fecal samples were collected from four patients (herein referred 
to as patient A, B, C and D) at the Department of Internal 
Medicine of the University Hospital Schleswig-Holstein, Campus 
Kiel, Germany (UK-SH). Patient A was treated with moxifloxacin 
(400 mg/day) for 13 days. Moxifloxacin is a fourth-generation 
synthetic fluoroquinolone antibacterial agent with a bactericidal 
effect inhibiting cell replication. AB treatment for patient B 
consisted of a combined therapy with penicillin G and clindamycin 
on the day of admission, and subsequendy with clindamycin alone 
(3x300 mg/day) for seven days. This semi-synthetic derivative 
belongs to the lincosamide class exerting a bacteriostatic effect due 
to the inhibition of protein synthesis. For patient C, AB therapy 
was initiated with cefazolin (3x2 g/ day) for seven days and 
continued with ampicillin/sulbactam (2x750 mg) for seven more 
days. Patient D received an amoxicillin (3 x 1000 mg/day) 
treatment. The antibiotics used for these two latter two patients 
belong to the beta-lactam class and have a bactericidal effect 
inhibiting cell envelope synthesis. Main features of patients and 
therapy are shown in Table 1. 

Fecal samples from patients (named A, B, C and D) were 
collected on the day of admission, before the antibiotic treatment 
(day 0), during and after AB therapy. In two cases (A and B), the 
last sample was taken 3 days after therapy, in the other two cases 
(B and D) the last sample was provided 7 and 28 days after 
treatment, respectively (Table 1). Patients did not present any 
intestinal disorder. Samples were collected in sterile tubes and 
stored at — 80°C until further processing. 

DNA extraction 

Tubes containing fecal samples with sterile PBS (containing, per 
liter, 8 g of NaCl, 0.2 g of KC1, 1.44 g of Na2HP04, and 0.24 g 
of KH2PQ4 [pH 7.2]) were centrifuged at 1250 g and 4°C for 2 



min to remove fecal waste. DNA was extracted from bacterial 
pellets using QIAamp® DNA Stool Kit (Quiagen) following the 
manufacturer's instructions. The product was concentrated by 
precipitation using 0.1 V of NaCl 3 M and 2 V of ethanol 100% 
and diluted in 75 (il of nuclease-free water. A standard agarose gel 
electrophoresis was run to check the integrity of DNA. The total 
DNA obtained was quantified with Nanodrop-1000 Spectropho- 
tometer (Thermo Scientific) and with the QuantiT PicoGreen 
dsDNA Assay Kit (Invitrogen). 

Amplification of the 16S rRNA gene 

For each sample a region of the 1 6S rRNA gene was amplified 
by polymerase chain reaction (PCR). The primers used were the 
universal E8F (5'-AGAGTTTGATCMTGGCTCAG-3') with 
adaptor A and 530R (5'-CCGCGGCKGCTGGCAC-3') with 
adaptor B using the sample-specific Multiplex Identifier (MID) for 
pyrosequencing. The amplified region comprises hyper-variable 
regions VI, V2 and V3. For each sample a 50 u.1 PCR mix was 
prepared containing 5 |il of Buffer Taq (10X) with 20 mM 
MgC12, 2 ul of dNTPs (10 mM), 1 ul of each primer (10 mM), 
0.4 |ll of Taq Fast start polymerase (5 u/ ul), 39.6 Jul of nuclease- 
free water and 1 n.1 of DNA template. PCR was run under the 
following conditions: 95° for 2 min followed by 25 cycles of 95° for 
30 s, 52° for 1 min and 72° for 1 min and a final extension step at 
72° for 10 min. The amplification process was checked by 
electrophoresis in agarose gel (1.4%). PCR products were purified 
using NucleoFast® 96 PCR Clean-Up Kit (Macherey-Nagel) and 
quantified with Nanodrop-1000 Spectrophotometer (Thermo 
Scientific) and with the QuantiT PicoGreen dsDNA Assay Kit 
(Invitrogen). The pooled PCR products were directiy pyrose- 
quenced. 

Total RNA extraction and double-strand cDNA synthesis 

Total RNA was extracted from fecal samples using RiboPur- 
e™-Bacteria Kit (Ambion). DNase treatment was applied to 
remove traces of genomic DNA from the eluted RNA using the 
same kit. The integrity of RNA was checked by electrophoresis 
in agarose gel (0.8%). The efficiency of the DNase treatment 
was checked by amplifying each RNA sample by PCR. To 
retro-transcribe total RNA into single-stranded cDNA the High- 
Capacity cDNA Reverse Transcription Kit (Ambion) was used. 
To synthesize double-stranded cDNA, 7.5 u.1 of Escherichia coli 
ligase buffer (10X), 2 ^1 of dNTPs (10 mM), 0.2 ul of E. coli 
RNAse H (5 u/ul), 3 ul of E. coli DNA pol I (10 u/ul), 0.5 ^1 of 
E. coli ligase (10 u/|il) and 41.8 \iX of nuclease-free water were 
added to each single-stranded cDNA sample. The mixture was 
placed in a Thermocycler at 15°C for 2 hours. Then, 2.5 |il of 
T4 DNA polymerase (3 u/(ll) were added and kept at 15°C for 
30 min. The metatranscriptome obtained thus was purified by 
precipitation and quantified using Nanodrop-1000 Spectropho- 
tometer (Thermo Scientific) and the QuantiT PicoGreen 
dsDNA Assay Kit (Invitrogen). A standard agarose gel electro- 
phoresis was run to check the integrity of double-stranded 
cDNA. 

Pyrosequencing 

For each sample, the total DNA (metagenome), double-stranded 
cDNA and amplicons of the 16S rRNA gene were sequenced with 
a Roche GS FLX sequencer and Titanium chemistry in the 
company Life Sequencing (V alencia, Spain) and in the Center for 
Public Health Research (CSISP-FISABIO) (Valencia, Spain). We 
obtained an average of 58,928, 41,838 and 4,872 reads per 
sample, respectively. 
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Table 1. Main features of the follow-up study. 



Patient Antibiotic 


Mode of action/ Antimicrobial effect 


Pathology 


Sampling date 


Samples 


A Moxifloxacin 


Cell replication inhibitor/ Bactericidal 


Bronchitis, pneumonia 


dayO-before AB 


A_before 








day3-during AB 


A3_D 








day6-during AB 


A6_D 








daylu-during AB 


A10_D 








day13-during AB 


A13_D 








3 days after AB 


A_after 


B Clindamycin 


Protein synthesis inhibitor/ Bacteriostatic 


Erysipelas 


dayO-before AB 


B_before 








day2-during AB 


B2_D 








day5-during AB 


B5_D 








day6-during AB 


B6_D 








28 days after AB 


B_after 


C Cefazolin/ Ampicillin/ Sulbactam 


Cell envelop synthesis inhibitor/ Bactericidal 


Infection pacemaker 


dayO-before AB 


C_before 








day3-during AB 


C3_D 








day6-during AB 


C6_D 








daylO-during AB 


C10_D 








3 days after AB 


C_after 


D Amoxicillin 


Cell envelop synthesis inhibitor/ Bactericidal 


Chronic sinusitis maxillans 


dayO-before AB 


D_before 








day3-during AB 


D3_D 








7 days after AB 


D_after 



AB, antibiotic; D, during the treatment. 
doi:1 0.1 371 /joumal.pone.0080201 .t001 



Taxonomic assignment of 16S rRNA amplicons 

We have used the Ribosomal Database Project (RDP) 
pyrosequencing pipeline [24] to trim off the MID and primers 
and to obtain the taxonomic classification. Sequences with a phred 
quality score less than 20 (Q20) and short length (<250pb) were 
discarded. We considered only annotations that were obtained 
with a bootstrap value greater than 0.8, leaving the assignation at 
the last-well identified level and consecutive levels as unclassified 
(uc). 

Taxonomic assignment of 16S rRNA transcripts 

Due to the procedure followed to obtain the metatranscriptome, 
the vast majority of transcripts belonged to ribosomal genes (16S 
and 23S). The 16S rRNA reads were obtained from the total 
cDNA by comparing the total reads against the Small Subunit 
rRNA Reference Database (SSUrdb) [25] with BLASTN [26] and 
an e-value of 10- . All sequences with detected homology were 
considered as 16S rRNAs and used to evaluate the phylogenetic 
diversity of the active bacteria. The taxonomic classification was 
performed in the same way as the amplicons. 

Analysis of total and active microbiota 

To study the phylogenetic structure of the bacterial community 
we applied two approaches that involved the 16S rRNA gene. The 
widely used analysis of 16S rRNA gene amplicons shows the 
composition of the total microbiota (16S rRNA gene). However, 
since the growing (active) bacteria contain more ribosomal RNA 
than latent or starved cells, studying the 16S ribosomal RNA 
transcripts enabled the active members of the microbiota to be 
identified (16S rRNA transcripts) [12,22]. We calculated sample 
diversity of the throughout the treatment for total and active 
bacteria by applying three parameters: two richness estimators, 
Chaol [27] and the abundance-based coverage (ACE) [28], and 



the Shannon index [29]. These estimators are implemented in 
package Vegan [30] under R software (http://cran.r-project.org) 
[31]. The biodiversity index and richness estimators were 
calculated after sub-sampling with the multiple_rarefactions.py 
script of QIIME to avoid the bias of the sequencing effort [32]. We 
used heat maps based on taxonomic composition to study the 
similarity between samples due to the relative abundance of each 
taxon using the Vegan library in the R software (http://cran.r- 
project.org) [30,31]. Canonical correspondence analysis (CCA) 
was performed to determine the relation between the sample 
composition and the class and mode of AB-treatment. To 
statistically assess the effect of such factors on the bacterial 
composition a multivariate ANOVA based on dissimilarity tests 
(Adonis) was applied, as implemented in the package Vegan, R 
software (http://cran.r-project.org) [30,31]. 

Metagenomics: functional analysis 

To eliminate reads that were artifact replicates of pyrosequenc- 
ing, we used the 454 Replicate Filter Program [33] with the 
following parameters: sequence identity cutoff = 1; length 
difference requirement = 0; number of beginning base pairs to 
check = 10. Unique reads were compared against the human 
genome using BLASTN [26] with an e-value of 10" 10 in order to 
remove human sequences. To identify the sequences encoding the 
ribosomal gene 16S rRNA we compared the dataset against the 
Small Subunit rRNA Reference Database (SSUrdb) described in 
Urich et al. [25] using BLASTN [26] with an e-value of 10" 1B . 
Sequences that did not give homology were used to identify the 
reads corresponding to the ribosomal gene 23S rRNA by 
BLASTN [26] against the Large Subunit rRNA Reference 
Database (LSUrdb) described in Urich et al. [25] with an e-value 
of 10~ 4 . Reads that matched with the LSUrdb were discarded. 
The remaining reads were compared to the NCBI-nr protein 
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database using BLASTX [26] to identify the protein-coding genes. 
Taxonomic assignment was based on the output of BLASTX 
applying the lowest common ancestor (LCA) algorithm. Fasta files 
were used to identify the Open Reading Frames (ORFs) by 
applying the facility of Fraggenscan from the web server of 
metagenomic analysis (WebMGA) [34] . To annotate the functions 
of the predicted ORFs, we applied HMMER 3.0 program [35] 
against TIGRFAM database [36] using default parameters. To 
identify the genes involved in resistance to antibiotics, we 
compared the identified ORFs against the Antibiotic Resistance 
Genes Database by BLASTp [37] with an e-value of 10" 10 . We 
used the ShotgunFunctionalizeR package [38] in the R software 
http://www.R-project.org/ [31] for functional comparison of 
metagenomes. Specifically, we applied the testGeneCategories.- 
dircomp function to compare the distribution of functional 
categories between groups of samples. The test is based on a 
Poisson model and compares each gene family of a higher 
functional category to decide if the category is statistically 
significant among two groups of samples [38] . 

Data accession number 

All sequences have been entered in the European Bioinfor- 
matics Institute database, under accession number ERP002192. 

Results 

Dynamics of total and active microbiota composition 
throughout therapy 

We analyzed total (16S rRNA gene) and active (16S rRNA 
transcripts) microbiota from the four patients (A, B, C and D) 
throughout AB treatment. The antibiotics administered to patients 
A, C and D had a bactericidal antimicrobial effect, whereas in 
patient B the effect was bacteriostatic. Regarding the mode of 
action, the antibiotic used in patient A was a cell replication 
inhibitor, in patient B it was an inhibitor of protein synthesis, 
whereas patients C and D were treated with a cell envelop 
synthesis inhibitor (Table 1). Each patient not only presented their 
own microbiota profile for both total (Figure 1A) and active 
(Figure IB) microbiota before treatment, but also there was 
apparently a rather large variation in bacterial taxa abundance 
during and after treatment, which we describe succinctly. In 
patient A, both total and active microbiota showed a high presence 
of the families Lachnospiraceae (Coprococcus and Roseburia genera) 
and Ruminococcaceae (Faecalibacterium, Blautia and Subdoligranulum 
genera) during AB treatment with fluoroquinolone (Figure 1). 
However, some genera such as Faecalibacterium and Subdoligranulum 
were negatively affected by the AB, while others such as Blautia, 
Coprococcus, Coprobacillus and Collinsella appeared to be resistant in 
the first stage of treatment. We also found that the bactericidal 
effect of AB had a negative impact on Bacteroides genus 
(Bacteroidetes phylum) in the first days of treatment, but the 
trend changed on day 13 with a great increase in its abundance. 
Treatment with clindamycin of patient B resulted in a high 
presence of Enterobacteriaceae [Escherichia, Salmonella genera), as 
shown in Figure IB. We also observed an increase in the Bacteroides 
genus after the 5th day of treatment in active microbiota. For 
patient C, Oscillibacteriaceae and Ruminococcaceae families 
(Firmicutes phylum) as well as Rikenellaceae and Bacteroidaceae 
(Bacteroidetes phylum) constituted the most abundant taxa in the 
total microbiota (Figure 1A). The first important change occurred 
on day 6 with an increase in Parabacteroides (Bacteroidetes phylum), 
which remained abundant after treatment. However, in the active 
microbiota we observed a shift towards the Bacteroidetes phylum 
(Alistipes and Bacteroides genera) at the two last time points. On the 



10th day of treatment, an increase in facultative anaerobic 
families, Enterobacteriaceae (Proteobacteria) and Enterococcaceae 
(Firmicutes) was found. Finally, in patient D, the initial microbiota 
composition consisted mainly of Enteriobacteriaceae (Escherichia 
genus) and Ruminococcaceae (Faecalibacterium genus). However, 
both genera were greatly affected by the antibiotic as there was an 
increase in resistant bacterial taxa of the Bacteroides genus 
(Bacteroidetes). 

After the AB course, patients A, C and D who received a 
bactericidal antimicrobial agent clustered together in both cases, 
total (Figure 2A) and active (Figure 2B) microbiota, apart from the 
patient treated with a bacteriostatic antibiotic (patient B) (Figure 1 
and 2). Moreover, we observed that the two patients treated with 
cell envelope synthesis inhibitors (C and D) grouped together in 
the case of the active microbiota (Figure 2B). 

For all patients, the diversity parameters (Shannon index, 
Chaol and ACE estimators) of both total and growing microbiota, 
showed notable fluctuations with a decrease in the number of 
bacterial taxa and evenness on the first days of treatment (Figure 
SI). At the end of the AB course, these three biodiversity 
parameters increased but they did not reach the initial values 
observed before AB therapy (Figure SI). 

Effect of the class of antibiotic 

To evaluate the pattern of variation shown by bacterial taxa or 
gene abundances and its relationship with two variables (the 
antimicrobial effect -bactericidal and bacteriostatic- and the mode 
of action of the antibiotic -protein synthesis inhibitor, cell 
replication inhibitor and cell envelope synthesis inhibitor-) we 
applied a CCA at the different levels: 16S rRNA gene, 16S rRNA 
transcripts, genes and the taxonomy of the identified coding 
regions (gene taxonomy). The results showed that these two factors 
(antimicrobial effect and mode of action) accounted for variability 
in a particular direction and with different strength (Figure 3). 
Figure 3 A shows that the first axis explained 19% of variability, 
splitting the total microbiota (16S rRNA gene) of the patients that 
were medicated with bactericidal AB (patients A, C and D) from 
the one treated with a bacteriostatic AB (patient B). The second 
axis explained 12% of variability, placing the two groups of 
samples treated with cell replication inhibitor (patient A) and 
protein synthesis inhibitor (patient B) antibiotics on one side of the 
graph; these ABs inhibit both essential and related cellular 
processes, such as DNA replication and protein synthesis. By 
contrast, the samples from patient C, treated with a cell envelope 
synthesis inhibitor AB affecting synthesis of the bacterial cell wall, 
fell on the other side of the graph. Both variables (antimicrobial 
effect and mode of action) introduced significant variance in the 
microbiota composition (Adonis test: p = 0.02, p = 0.04, 
respectively). 

Regarding active microbiota (Figure 3B), the first and second 
axes explained 12% and 6% of the total variability, respectively. 
With respect to the first axis, the samples from patient B (protein 
synthesis inhibitor antibiotic) situated on the right side of the 
graph. The latter AB, as occurred for the total microbiota, 
introduced higher variance. The second axis separated the samples 
taken from patients C and D, treated with cell envelope synthesis 
inhibitor AB from the rest. Despite both factors, the antimicrobial 
effect and mode of action were not significant (Adonis test: p = 
0.18, p= 0.069), the second explained more variability (Adonis 
test, p = 0.069). 

The CCA applied at gene level (Figure 3C) showed a 
distribution similar to that found for the total microbiota (Figure 
3 A), with the first and the second axes explaining 25% and 6% of 
the total variance, respectively. The samples are separated chiefly 
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Figure 1. Microbiota composition of patients A, B, C, and D. (A) total microbiota (16S rRNA gene) (B) active microbiota (16S rRNA transcripts). 

The mode of action for each AB used is indicated. 

doi:10.1371/journal.pone.0080201.g001 



by antimicrobial effect and then by mode of action. However, in 
this case the strength of the different vectors is weaker, probably 
due to the great functional homogeneity of the gut microbial 
community (Adonis test: antimicrobial effect p= 0.27, mode of 
action p = 0.41). 

Finally, we performed a CCA using the taxonomy of the 
identified coding regions (Figure 3D). The two axes explained 19% 
and 7% of the total variation of the data. As can be seen, the 
different classes of antibiotics affected the taxonomy of the 
identified coding regions in a similar way to the results reported 
for total microbiota (Figure 3A) and gene analysis (Figure 3C) 
(Adonis test: antimicrobial effect p = 0.044, mode of action p = 
0.049). 



Functional analysis of the metagenomes 

The functional annotation of the ORFs was derived using the 
TIGRFAM database, providing a hierarchical order: main roles, 
the highest functional level (described in Figure S2), sub-roles, 
more specific metabolic functions for each one of the main roles 
and genes [36] . Regarding the main roles, the profiles were fairly 
homogeneous for all patients and time points. The most abundant 
categories before, during and after treatment were "Protein 
synthesis", "Energy metabolism", "Cellular processes" and 
"Transport and Binding Proteins" with average values of relative 
abundance of 13.5%, 13.2%, 9.6% and 9.5%, respectively, which 
highlights the importance of these functions performed by the gut 
microbiota (Figure S2). However, when considering the sub-roles, 
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Figure 2. Heat map and clustering based on taxon composition and abundance. (A) total microbiota, (B) active microbiota. Colors in the 
figure depict the percentage range of sequences assigned to main taxa (abundance >1% in at least one sample). 
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we detected significant changes in the corresponding profiles for 
each patient before, during and after treatment. For all patients we 
detected a total of 53 sub-roles that differed significantly in gene 
content (Table 2). Only two sub-role functional categories changed 
significantly in all patients: "Menaquinone and ubiquinone" 
(within the main function "Biosynthesis of cofactors, prosthetic 
groups, and carriers") and "Carbohydrates, organic alcohols, and 
acids" (within the main function of "Transport and binding 
proteins"). Genes participating in the biosynthesis of menaquinone 
and ubiquinone were under-represented during the treatment for 
all patients, except in the case of patient B. Regarding the sub-role 
of "Carbohydrates, organic alcohols, and acids" the gene 
functions were overrepresented during the treatment in patients 
A, B and C and under-represented in patient D. Most of the genes 
belonging to this functional group were related to the phospho- 
transferase system (PTS), which is essential for translocating 
carbohydrates in bacteria [39] . Within this family, we have found 
genes involved in the transport of various sugars such as mannose, 
fructose, sorbose, glucitol or glucose. It is noteworthy that the 
related sub-role "PTS" (within the main role "Signal transduc- 
tion"), associated to genes participating in regulation, was also 
over-represented during the treatment in patients A and C and 
under-represented in patient D. 

The changes in the sub-roles "Biosynthesis and degradation of 
surface polysaccharides and lipopolysaccharides" and "Other" 
(from the "Cell envelope main role") were significant in patients A, 
B and D. In general, we detected a lower presence during 
treatment of genes involved in the synthesis of lipopolysaccharides 
(LPS). The sub-role "Pathogenesis" (within the main category 
"Cellular processes") and "Degradation" (within the "Fatty acid 
and phospholipid metabolism" main category) decreased signifi- 
cantly during treatment for patients A and D and increased in 
patient B. The fatty acid and phospholipid metabolism genes were 
involved in fatty acid beta-oxidation. On the contrary, the sub-role 



"Sporulation and germination" (within "Cellular processes") was 
more abundant during treatment in patients A and D, with most of 
the genes being involved in different stages of endospore 
formation. Finally, we found that the sub-role "TCA cycle" 
(within the main category "Energy metabolism") and "Amino 
acids, peptides and amines" (within "Transport and binding 
proteins) underwent an increase in the number of genes encoding 
different enzymes of the citric acid cycle during antibiotic 
treatment for patient A and B and a decrease for D. Regarding 
the transport of amino acids, peptides and amines, we found the 
presence of genes encoding ABC transporters for amino acids and 
urea. 

Analysis of the resistome 

We performed a search of the resistome by identifying AB 
resistance genes in the 19 metagenomes analyzed by comparing 
the predicted ORFs against the Antibiotic Resistance Genes 
Database [37]. We identified the resistance genes that represented 
0.2%, 0.8%, 0.22%, and 0.5% of the total determinants found for 
patients A, B, C and D respectively. We found that while patients 
A, B, and C showed an increase in resistance genes after 
treatment, patient D presented the lowest relative abundance of 
these determinants, decreasing from 0.81% to 0.14% (Figure 4A). 
Figure 4B shows the profiles of resistance genes that varied during 
the treatment for patient A, B and C, administered with antibiotics 
belonging to different classes: fluoroquinolone, lincosamide and 
beta-lactams, respectively (Table 1). Overall, we observed that the 
resistance induced by each antimicrobial was associated with other 
resistance determinants. Also, we found that its profiles matched 
fairly well with the oscillatory dynamics of the surviving bacterial 
community. Patient A showed an increase in the relative 
abundance of the total resistance genes at the end of the 
treatment, raising values from 0.18% before treatment to 0.28% 
after the AB course. Fluoroquinolone resistance, multidrug 
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Figure 3. Canonical Correspondence Analysis (CCA) of the patients A, B, C and D in the follow-up study. (A) total microbiota, (B) active 
microbiota, (C) genes and (D) gene taxonomy. The antimicrobial effect is represented as a vector with two levels (bactericidal and bacteriostatic). The 
mode of AB action is represented as a vector with three levels (cell envelop synthesis inhibitor, cell replication inhibitor and protein synthesis 
inhibitor). 
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resistance efflux pump, appeared on day 10, when the microbiota 
composition was dominated by the Firmicutes phylum, with high 
abundance of members of the genus Entemcoccus, described as 
resistant to this type of antimicrobial [40] . This patient presented 
high relative abundance of tetq gene, which confers tetracycline 
resistance before and after the treatment, being Bacteroides genus 
one of the most abundant taxa. Bacitracin profile showed a 
maximum on day 6, with baca gene being associated to Streptococcus 
and Clostridiales taxa (Figure 1 and 2). Patient B showed a strong 
increase in the relative abundance of resistance genes during 
treatment, increasing from 0.29% up to 0.89%. In fact all the 



genes increased in abundance after AB treatment except those 
involved in bacitracin resistance. The most remarkable increase 
was found in a group of genes coding for multidrug resistance 
efflux pump, which confer resistance against clindamycin and 
related antimicrobials (aminoglycoside, glycylcycline, beta-lactam, 
macrolide, and acriflavine). Patient C also showed an increase in 
the total relative abundance of resistance genes at the last time 
point and after treatment, from 0.16% to 0.36%. The genes that 
code for multidrug resistance efflux pumps are the most abundant 
on the 10 th day of treatment. However, beta-lactamase genes 
increased throughout the AB course and reached the maximum at 
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Table 2. Functional profiles. 



Patient 


Main Role 


Sub-Role 


A 


B 


c 


D 


Amino acid biosynthesis 


Glutamate family 


f 2.06E-004 


NS 


NS 


NS 




Histidine family 


f 0.02 


NS 


NS 


NS 




Serine family 


f 4.36E-004 


NS 


NS 


f 0.01 


Biosynthesis of cofactors* 


Biotin 


NS 


t 0.02 


NS 


I 2.33E-004 




Glutathione and analogs 


NS 


1 2.95E-003 


NS 


I 0.04 




Menaquinone and ubiquinone 


I 0.02 


f 2.19E-005 


I 0.04 


| 4.68E-003 




Molybdopterin 


NS 


NS 


NS 


I 4.35E-003 




Pantothenate and coenzyme A 


T 0.01 


NS 


NS 


NS 




Other 


f 1.33E-004 


NS 


NS 


NS 


Cell envelope 


Biosynthesis and degradation of surface** 


I 0.05 


f 3.01E-005 


NS 


4 0.02 




Other 


f 0.01 


I 0.01 


NS 


f 0.02 




Surface structures 


NS 


NS 


NS 


T 0.04 


Cellular processes 


Biosynthesis of natural products 


NS 


NS 


NS 


I 0.05 




Cell division 


NS 


I 0.01 


NS 


1 3.07E-004 




Chemotaxis and motility 


NS 


NS 


NS 


T 0.01 




Detoxification 


NS 


NS 


NS 


1 0.04 




DNA transformation 


I 2.32E-003 


NS 


I 0.01 


NS 




Pathogenesis 


I 9.43E-004 


T 1.47E-005 


NS 


I 2.22E-009 




Sporulation and germination 


f 0.03 


I 7.07E-017 


NS 


1 2.65E-017 




Toxin production and resistance 


NS 


I 2.85E-004 


NS 


NS 


Central intermediary metabolism 


Amino sugars 


NS 


NS 


NS 


1 0.03 




Nitrogen metabolism 


NS 


t 0.02 


NS 


NS 


DNA metabolism 


Chromosome-associated proteins 


NS 


NS 


J, 0.01 


NS 




Restriction/modification 


f 0.02 


NS 


NS 


NS 


Energy metabolism 


Aerobic 


I 0.01 


NS 


NS 


J, 0.01 




Amino acids and amines 


NS 


f 9.26E-007 


NS 


I 2.22E-009 




Anaerobic 


NS 


NS 


NS 


I 4.80E-005 




Biosynthesis and degradation of polysaccharides 


I 1.58E-003 


NS 


NS 


NS 




Chemoautotrophy 


NS 


1 0.01 


NS 


NS 




Electron transport 


NS 


t 0.01 


NS 


NS 




Entner-Doudoroff 


NS 


T 2.85E-004 


NS 


1 0.01 




Fermentation 


NS 


NS 


NS 


T o.oi 




Pentose phosphate pathway 


NS 


T 0.01 


NS 


1 0.03 




Sugars 


NS 


T 2.47E-005 


NS 


1 0.04 




TCA cycle 


f 0.05 


T 0.02 


NS 


I 4.59E-003 


Fatty acid and phospholipid metabolism 


Degradation 


I 0.05 


| 0.04 


NS 


| 0.03 


Protein fate 


Protein and peptide secretion and trafficking 


NS 


NS 


NS 


I 0.02 




Protein folding and stabilization 


NS 


I 0.01 


NS 


NS 


Protein synthesis 


Other 


I 0.03 


I 4.57E-005 


NS 


NS 




Ribosomal proteins: synthesis and modification 


NS 


I 2.85E-004 


NS 


f 4.80E-005 




tRNA and rRNA base modification 


NS 


NS 


t 2.04E-003 


NS 


Regulatory functions 


Other 


NS 


NS 


t 0.01 


NS 


Signal transduction 


PTS 


t 1.40E-005 


NS 


t 3.59E-009 


I 7.09E-009 




Two-component systems 


NS 


NS 


NS 


I 1 .52E-003 


Transcription 


DNA-dependent RNA polymerase 


NS 


1 4.49E-003 


NS 


t 0.03 


Transport and binding proteins 


Amino acids, peptides and amines 


t 4.36E-004 


t 2.85E-004 


NS 


I 4.32E-005 




Anions 


I 4.36E-004 


NS 


NS 


NS 




Carbohydrates, organic alcohols, and acids 


t 4.36E-004 


t 3.62E-006 


t 1.46E-004 


I 3.76E-033 
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Patient 


Main Role 


Sub-Role 


A 


B 


C 


D 




Cations and iron carrying compounds 


t 3.03E-003 


NS 


NS 


1 0.04 




Nucleosides, purines and pyrimidines 


NS 


t 2.61E-003 


NS 


J, 1.52E-003 




Other 


NS 


t 1.33E-003 


NS 


I 5.56E-004 




Porins 


NS 


t 0.01 


NS 


NS 


Unknown function 


Enzymes of unknown specificity 


I 0.04 


NS 


NS 


I 4.97E-005 



Main roles and sub-roles that change significantly during treatment and their associated p-values (p-value < 0.05). The upward arrow indicates those categories that 
were more abundant during treatment and the downward arrow those that were less abundant. NS, not significant. 
*Biosynthesis of cofactors, prosthetic groups, and carriers. 

**Biosynthesis and degradation of surface polysaccharides and lipopolysaccharides. 
doi:1 0.1 371 /journal.pone.0080201 .t002 



the end of treatment. Patient C also presented a high abundance of 
tetracycline resistance genes before treatment {teto, tetq and tetw) 
associated to different taxa (Blautia, Bactemides, Clostridium, and 
Ruminococcus), which have been described as resistant to this 
antibiotic [41,42,43] but underwent a dramatic decrease on day 
10 associated with a major presence of the Proteobacteria phylum. 

Discussion 

Dynamics of the gut microbiota structure over the AB 
course 

The human gut microbiota consists of a highly complex 
community whose members establish close relationships with the 
host. ABs have strong direct and indirect effects on the human gut 
microbiota and consequendy on the functions they perform, 
affecting the ecosystem maintenance and therefore host physiology 
[7,13]. The microorganisms that carry certain genetic determi- 
nants have an advantage under AB pressure, allowing them to 
survive and grow. It is well known that the human gut microbiota 
presents a high inter-individual variability and that its composition 
depends on factors such as genetics, age, diet, health status and 
AB-therapy, among others. 

In our study, each patient presented their own initial microbiota 
and thus there was an individual response to AB treatment with 
fluctuations in the bacterial diversity and composition for both 
total and active gut microbiota. These results highlight the 
importance of the initial microbial structure in shaping the 
changes in microbiota during the AB course. The individual 
character of the response and incomplete recovery of initial 
microbiota after AB treatment has previously been described by 
Dethlefsen and coworkers [18] in a follow-up study of three 
patients that received two courses of ciprofloxacin. However, we 
also observed that the selection of resistant microorganisms led to a 
similar microbiota composition after analogous antibiotic treat- 
ment. Thus, AB seems to have a major impact on the structure of 
the final bacterial community. 

The gut microbiota has been described as an ecosystem that is 
relatively resistant to perturbations [16]. However we observed 
that a particular assembly of microorganisms can confer greater 
resistance to a disturbance than others in terms of presence and 
abundance of taxa, which could be related with the specific effect 
of the AB. In patient A, there was a decrease in Faecalibacterium and 
Bactemides genera during the AB course, with AB-resistant strains 
appearing at the end of treatment. However, on the first days of 
treatment, other butyrate-producing taxa (Roseburia and Lachnospir- 
aceae incertae sedis) and H2-consuming bacteria [Blautia, Collinsella 



and Bifidobacterium) were present as active microbiota, obtaining 
energy sources for the colonocytes of the host. A similar behavioral 
pattern of these members of the intestinal microbiota has also been 
reported by the above mentioned group of Dethlefsen et al. 
[15,18] when they used ciprofloxacin, as in our patient A, an AB 
belonging to the cell replication inhibitor group. In the case of 
patient B, since clindamycin affected anaerobic bacteria, there 
were marked decreases in Bactemides and Blautia genera in the 
active microbiota just after AB administration (Figure 2B). 
However, three days later, high abundances of Bactemides were 
detected in the active bacterial community, suggesting these 
bacteria acquired resistance. The presence of clindamycin- 
resistant Bactemides in gut microbiota has also been described in 
other studies [9,44]. Moreover, we observed that the decrease in 
anaerobic bacteria is compensated for by an important increase in 
members of the family Enterobacteriaceae. Patients C and D were 
treated with ABs which have a similar mode of action as both are 
of the (3-lactam class. As stated, the initial microbiota composition 
was very different between both subjects, showing a differential 
response to ABs. However, the active microbiota changed 
throughout AB treatment with both patients acquiring a similar 
composition by the end of it (Figure 2B). Patient C received a 
combination of two ABs, Cefazolin and Ampicillin/ Sulbactam, 
which cover a broad spectrum of microorganisms and showed a 
significant increase in Parabacteroides and Bactemides genera. 
Interestingly, resistance genes against ampicillin and cephalosporin 
in these two taxa have been described previously [43]. On day 10 
of treatment, an increase in the Enterobacteriaceae family 
occurred and some of its genera, such as Escherichia or Klebsiella 
are considered as opportunistic pathogens [45], suggesting that AB 
use creates opportunistic infections by these harmful microorgan- 
isms. Patient D was treated with amoxicillin, described as active 
against some Proteobacteria such as Escherichia or Klebsiella. During 
treatment, these genera were almost eliminated, whereas Bactemi- 
des, Blautia and Faecalibaterium taxa proved less susceptible to 
treatment, as occurred in patient C. It is worth pointing out that 
this bacterial profile, with Bactemides, Blautia and Faecalibaterium 
after the AB stress, has been found in the case of bactericidal 
agents but not when a bacteriostatic antimicrobial was used. 
Furthermore, a previous study [12] found a similar pattern of 
active bacteria after beta-lactam treatment. 

Antimicrobial effect and mode of action of ABs on the 
gut microbiota 

It has been stated that external variables such as ABs shift the 
microbial composition [46]. In our study, the class of AB 
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Figure 4. Resistance gene profiles. (A) The dashed lines represent the relative abundance of the total number of resistance genes for patients A, 
B, C, and D. (B) Relative abundance of the resistance genes throughout AB treatment for patients A, B, and C. The symbol "*" highlights the resistance 
gene profiles which coincide with the antibiotic administered to patients C, A and B, respectively. 
doi:1 0.1 371/journal.pone.0080201 .g004 



significantly shaped the microbiota on the basis of the antimicro- 
bial effect (bactericidal or bacteriostatic) and the mode of action. 
In addition, we have found that specific mechanisms of action 
affect some organisms more than others, leading the bacterial 
community towards an alternative temporary equilibrium state. 
Clindamycin (protein synthesis inhibitor) introduced higher 
variance in microbiota composition than the other agents, giving 
way to a different bacterial community structure. This was 
probably due to the bacteriostatic nature of clindamycin when 
compared to the bactericidal effect of the other AB treatments. 
Interestingly, in the case of the active microbiota representing the 
surviving community, the bacterial composition was affected by 
the mode of action rather than the antimicrobial effect clearly 
distinguishing the three modes of action (Figure 3B). At a 
functional level, the microbial community profile was driven by 
the antimicrobial effect rather than by the mode of action. 
However, the strength of the AB class, considered as an external 
factor exerted at gene level was less intense, resulting in major 
uniformity. 

AB impact on bacterial metabolic functions 

High homogeneity was observed in the main roles for all the 
patients. This uniformity at a functional level has been also shown 
in both DNA and RNA-based surveys [22,47,48,49,50] since the 
microbiota is characterized by high functional redundancy. When 
we considered sub-roles, 51% with significant variation corre- 
sponded to inter-individual variability, representing the specific- 
subject response to AB course. The over-representation of genes 
involved in sugar transport in most of the patients suggested that 
this functional category could play an important role under stress 
conditions, as is the case of AB treatment. The phosphotransferase 
system (PTS), in addition to its main role in sugar transport, which 
is an essential function in itself, is involved in different regulatory 
processes such as stress response in bacteria and, hence, it could 
confer some extra advantages in presence of ABs [51]. Then, an 
efficient system of importing sugars could facilitate the energetic 
metabolism and, therefore, it could counteract the negative effect 
of ABs on the bacterial growth. Perez-Cobas and coworkers [12] 
showed an increase in proteins belonging to the glycolysis pathway 
and pyruvate metabolism, as well as higher expression of genes 
related to energy metabolism/ sugars category during beta-lactam 
treatment. 

As we mentioned previously, the bacteriostatic effect drives 
the bacterial community to a characteristic composition, which 
is also reflected at a functional level. In patient B, most of the 
functional categories over-represented during treatment could 
be related with the increase in Enterobacteriaceae members. In 
this regard, we found an increase in the number of genes 
involved in lipopolysaccharide synthesis, which is the main 
component of the outer membrane for most Gram-negative 
bacteria. This barrier plays an important role in nutrient 
uptake and also confers resistance against ABs [52]. Likewise, 
the genes of secretion systems typical of Gram negative 
bacteria pathogenesis showed an increase only in patient B. 
However, patients A, C, and D, who received bactericidal 
treatment, presented a high abundance of Gram positive 
Ruminococcaceae and Lachnospiraceae families, and an over- 
representation of genes involved in endospore formation, a 
resistance mechanism typical of Gram positive bacteria. 



Another category that presented differences between bacteri- 
cidal and bacteriostatic ABs was catabolism of fatty acids and 
phospholipids, to produce acetyl-CoA through the beta- 
oxidation process. As clindamycin inhibits mainly the anaer- 
obic bacteria, the genes belonging to this sub-role were more 
abundant in patient B whose bacterial composition proved rich 
in Enterobacteriaceae family members. 

Changes in the resistome 

It has previously been pointed out that the AB usage is the 
most influential agent in the spread and stabilization of 
resistance genes in the gut environment [13]. One of the 
multidrug resistant genes that increased in patients A and B 
was a multidrug resistance efflux pump, which confers 
resistance against aminoglycoside, glycylcycline, beta-lactam, 
macrolide and acriflavine antibiotics. Since these ABs have 
different properties such as spectrum or mode of action, the 
transmission of these genes to a pathogen could hinder clinical 
treatment in the event of infection. These resistance genes have 
been described in some Proteobacteria genera such as 
Escherichia or Klebsiella, thereby supporting the increase in the 
abundance of these genera during treatment, principally in 
patient B. Patients A, B and C reached higher values of gene 
resistance abundances after AB treatment, with patient B, who 
was treated with clindamycin, attaining maximum values. It 
has been reported that besides the strong effect on the 
microbial composition, clindamycin also promotes increased 
AB resistance, which can persist in the microbial population for 
a long time [53]. In contrast, patient D showed a decrease in 
the relative abundance of resistance genes in the bacterial 
community. In fact, we found different dynamics in patient D 
as compared to the other three patients. This sample presented 
an initial composition with prevalence of the Enterobacteria- 
ceae family, which has been described as a considerable source 
of resistance genes [54] and hence, the data indicated that 
these taxa were strongly affected by the ABs. Thus, the final 
resistome in the human gut after AB therapy would be 
determined by the resistance genes carried by the surviving 
bacteria and by the class of AB administered. 

Conclusions 

In this study, using high-throughput methodology, we have 
provided new insights into the complex antibiotic resistance 
scenario, related to the different modes of action of antibiotics and 
the consequences for the gut microbiota composition and function 
during antibiotic therapy. We have shown that specific properties 
of ABs such as antimicrobial effects or mode of action, are 
powerful forces for the selection of intestinal microbiota, and are 
partially responsible for the shifts in bacterial composition during 
AB therapy. The resulting structure of the microbial community 
showed its specific metabolic capabilities giving a different 
functional profile. Additionally, we have shown that the AB also 
modified the resistome composition, increasing the abundance of 
resistance genes in the gut environment, which is also important in 
shaping the post-treatment composition of the microbiota. 
However, further research into a larger group of subjects would 
be necessary to establish a quantitative evaluation of changes in 
gut microbiota. 
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